      subroutine ornl_dfir2d(    f,incf,ldf,ifx0,n_fx,ify0,n_fy,
     $			    g,incg,ldg,igx0,n_gx,igy0,n_gy,
     $			    h,inch,ldh,ihx0,n_hx,ihy0,n_hy,
     $			    alpha, beta)
      double precision f(ldf,*)
      double precision g(ldg,*)
      double precision h(ldh,*)
      double precision alpha, beta
      integer incf, ldf, ifx0, n_fx, ify0, n_fy
      integer incg, ldg, igx0, n_gx, igy0, n_gy
      integer inch, ldh, ihx0, n_hx, ihy0, n_hy
c-----------------------------------------------------------------------------
c
c   dfir2d  performs a 2D convolution in the time domain :
c	h(i,j) = beta * h(i,j) + alpha * Sum.Sum[ f(k,l) * g(i-k,j-l) ]
c
c-----------------------------------------------------------------------------
c
c   PARAMETERS:
c
c	f:	Vector containing sequence "f"
c	incf:	Increment between two successive lines   of "f"
c	ldf:	Increment between two successive columns of "f"
c	ifx0:	Index of the first element of each column of "f"
c	n_fx:	Number of elements (lines) of each column of "f"
c	ify0:	Index of the first column of "f"
c	n_fy:	Index of the  last column of "f"
c	
c	g:	Vector containing sequence "g"
c	incg:	Increment between two successive lines   of "g"
c	ldg:	Increment between two successive columns of "g"
c	igx0:	Index of the first element of each column of "g"
c	n_gx:	Number of elements (lines) of each column of "g"
c	igy0:	Index of the first column of "g"
c	n_gy:	Index of the  last column of "g"
c	
c	h:	Vector containing sequence "h"
c	inch:	Increment between two successive lines   of "h"
c	ldh:	Increment between two successive columns of "h"
c	ihx0:	Index of the first element of each column of "h"
c	n_hx:	Number of elements (lines) of each column of "h"
c	ihy0:	Index of the first column of "h"
c	n_hy:	Index of the  last column of "h"
c
c	alpha: Scaling factor for the Convolution
c	beta:  Scaling Factor for the Output on Entry
c
c-----------------------------------------------------------------------------
c
c   PLEASE NOTE:
c	Please note that the array pointers must all point to the first 
c	element of the array "(ifx0,ify0)", "(igx0,igy0)" and "(ihx0,ihy0)"
c	 If "f" for example is defined as
c	    dimension f(-25:45,10:21)
c	Then "dfir2d" must be called with the following parameters
c	    call dfir2d( f(-25,10),(45-(-25)+1),-25,45,10,21 ... )
c
c-----------------------------------------------------------------------------
c
c   This Fortran Subroutine written by
c	Jean-Pierre Panziera
c	Silicon Graphics Inc
c	September 27, 1991
c
c-----------------------------------------------------------------------------

      double precision zero, one
      parameter (zero = 0.0, one = 1.0)
      integer iy0, iy1, ify1, igy1, ihy1
      integer k1, k2, k, j, ix, iy, jf, jg, jh

c-----------------------------------------------------------------------------
c
c   Compute Y Boundaries
c
c-----------------------------------------------------------------------------
	ify1 = ify0 + n_fy - 1
	igy1 = igy0 + n_gy - 1
	ihy1 = ihy0 + n_hy - 1

	iy0 = ihy0
	if( iy0 .lt. (ify0+igy0))	iy0 = ify0+igy0
	iy1 = ihy1
	if( iy1 .gt. (ify1+igy1))	iy1 = ify1+igy1

c-----------------------------------------------------------------------------
c
c   Scale the output
c
c-----------------------------------------------------------------------------
      if( beta .eq. zero) then
	do j = 1, n_hy
	    ix = 1
	    do k = 1, n_hx
		h(ix,j) = zero
		ix = ix + inch
	    end do
	end do
      else if( beta .ne. 1.0) then
	do j = 1, n_hy
	    ix = 1
	    do k = 1, n_hx
		h(ix,j) = beta * h(ix,j)
		ix = ix + inch
	    end do
	end do
      endif

      if( alpha .eq. zero ) return

c-----------------------------------------------------------------------------
c
c   Call the 1D convolution 
c
c-----------------------------------------------------------------------------
	do j = iy0, iy1
	    k1 = max( igy0, j-ify1)
	    k2 = min( j-ify0, igy1)
	    jh = j-ihy0+1
	    do k = k1, k2
		jf = j-k-ify0+1
		jg = k-igy0+1
		call dfir1d(	f(1,jf), incf, ifx0, n_fx,
     $				g(1,jg), incg, igx0, n_gx,
     $				h(1,jh), inch, ihx0, n_hx,
     $				alpha, one)
	    end do
	end do
c-----------------------------------------------------------------------------
      return
      end


